MODULE inout
  !-----------------------------------------------------------
  ! File utilities for the KRUUNUNHAKA basis set tool kit
  !-----------------------------------------------------------
  ! See CHANGELOG
  !-----------------------------------------------------------
  ! Copyright (C) 2003-2008 Pekka Manninen, 2010 Jussi Lehtola
  !
  ! This program is distributed under the terms of the GNU General
  ! Public License as published by the Free Software Foundation;
  ! either version 2 of the License, or (at your option) any later version.
  !
  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
  !
  ! Please note that this program as well as its 
  ! author must be properly cited whenever the program or some parts 
  ! originated on it are employed. 
  !
  !-----------------------------------------------------------

  USE error
  USE textutil
  USE definitions
  IMPLICIT NONE
CONTAINS
 !-----------------------------------------------------------
  ! Read in the number of exponents and contractions in Dalton format
  SUBROUTINE read_dalton_nbas(filename, nbas)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    CHARACTER(LEN=MAXLEN) :: line
    CHARACTER(LEN=MAXLEN) :: dummy
    INTEGER, DIMENSION(MAXL,2), INTENT(OUT) :: nbas
    INTEGER :: sta, i

    OPEN(10, FILE=filename, STATUS='OLD', ACTION='READ', IOSTAT=sta)
    IF (sta /= 0) THEN
       WRITE (*,*) 'Error opening file ', filename,&
            ' , error code:', sta
       RETURN
    END IF

    nbas = 0
    nbasloop: DO 
       ! Read in line
       READ(10,'(A)', IOSTAT=sta) line
       IF (sta /= 0) THEN
          EXIT nbasloop 
       END IF

       ! Look for angular momentum shells
       DO i = 1, MAXL
          IF (str_lo_to_up(line(3:3)) == l_name(i)) THEN
             ! Shell found, read in the number of contractions
             READ(10,*) dummy, nbas(i,1), nbas(i,2)
          END IF
       END DO
    END DO nbasloop
    CLOSE(10)

    IF(MAXVAL(nbas)==0) THEN
       CALL error_stop('Unable to find any basis functions in Dalton format!')
    END IF

  END SUBROUTINE read_dalton_nbas
 !-----------------------------------------------------------
  ! Read in basis set in Dalton format
  SUBROUTINE read_dalton_coeffs(filename, nbas, X_l, C_l)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    CHARACTER(LEN=MAXLEN) :: line
    INTEGER, DIMENSION(MAXL,2), INTENT(IN) :: nbas
    INTEGER :: sta, i, j, k
    REAL(KIND=prec), DIMENSION(:,:), INTENT(OUT) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(OUT) :: C_l

    OPEN(10, FILE=filename, STATUS='OLD', IOSTAT=sta, ACTION='READ')
    IF (sta /= 0) THEN
       WRITE (*,*) 'Error opening file ', filename,&
            ' , error code:', sta
       RETURN
    END IF

    ! Zero out arrays
    X_l=0.0_prec
    C_l=0.0_prec

    readxloop: DO
       ! Read the next line in the input file
       READ(10,'(A)', IOSTAT=sta) line
       IF (sta /= 0) THEN
          EXIT readxloop 
       END IF

       !! Look for new angular momentum shells
       ! Loop over l values
       DO i = 1, MAXL
          ! If the input line contains a shell, then
          IF (str_lo_to_up(line(3:3)) == l_name(i)) THEN
             ! Skip the next line, which contains the numbers of contractions
             READ(10,*) line
             DO j = 1, nbas(i,1)
                ! Read in the exponents and the coefficients
                READ(10,'(A1500)') line
                READ(line,*) X_l(i,j), (C_l(i,j,k), k=1,nbas(i,2))
             END DO
          END IF
       END DO
    END DO readxloop
    CLOSE(10)
    RETURN
  END SUBROUTINE read_dalton_coeffs
 !-----------------------------------------------------------
  ! Write basis in Dalton format
  SUBROUTINE write_dalton_basis(filename,nbas,X_l,C_l)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    CHARACTER(LEN=20) :: printform
    REAL(KIND=prec), DIMENSION(:,:) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:) :: C_l
    INTEGER, DIMENSION(MAXL,2), INTENT(IN) :: nbas
    INTEGER :: l, j, k

    IF(SIZE(C_l,2) /= SIZE(X_l,2)) THEN
       WRITE(*,*) 'Error in write_dalton_basis: contraction', &
            ' coefficients do not correspond to exponents!'
       CALL abort()
    END IF
    OPEN (20, FILE=filename, STATUS='replace')
    DO l = 1, MAXL
       WRITE(printform,'(A9,I0,A6)') '(G16.8,A,', nbas(l,2), 'G16.8)'
       IF(nbas(l,1) /= 0) THEN
          WRITE(20,'(A,A,A)') '! ', l_name(l), &
               '-TYPE FUNCTIONS'
          WRITE(20,'(A,I4,I4)') 'H', nbas(l,1), nbas(l,2)

          DO j = 1, nbas(l,1)
             WRITE(20,printform) X_l(l,j), ' ', &
                  (C_l(l,j,k), k=1,nbas(l,2))
          END DO
       END IF
    END DO
    CLOSE (20)      
  END SUBROUTINE write_dalton_basis
 !-----------------------------------------------------------
  ! Find out number of exponents and coefficients from basis in 
  ! Gaussian'94 format
  SUBROUTINE read_gaussian_nbas(filename,nbas)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    CHARACTER(LEN=MAXLEN) :: line
    CHARACTER(LEN=MAXLEN) :: dummystr
    INTEGER, DIMENSION(MAXL,2), INTENT(OUT) :: nbas
    INTEGER :: sta
    INTEGER :: i, j, l
    INTEGER :: nxtmp
    INTEGER, DIMENSION(MAXL) :: nxmax
    ! Temporary arrays to place the exponents and coefficients in
    REAL(KIND=prec), DIMENSION(:,:), ALLOCATABLE :: exps
    REAL(KIND=prec) :: expn, dummyreal
    LOGICAL :: found

    ! Initialize nbas
    nbas=0

    ! First, find out an upper limit for the number of exponents
    nxmax=0
    OPEN (10, FILE=filename, STATUS='old', ACTION='READ', IOSTAT=sta) ! Open input file

    IF (sta /= 0) THEN
       WRITE (*,*) 'Error opening file ', filename,&
            ' , error code:', sta
       RETURN
    END IF

    READ (10,*) line ! Skip first line

    nxloop: DO 
       READ(10,'(A)', IOSTAT=sta) line
       IF (sta /= 0) THEN
          EXIT nxloop
       END IF

       DO l = 1, MAXL ! Look for angular momentum shells
          IF (str_lo_to_up(line(1:1)) == l_name(l)) THEN ! Found one
             ! Read in the number of exponents
             READ(line,*) dummystr, nxtmp, dummyreal
             ! Increment nxmax
             nxmax(l) = nxmax(l) + nxtmp
             ! Skip the next nxtmp lines
             DO i=1,nxtmp
                READ(10,'(A)') dummystr
             END DO
          END IF
       END DO
    END DO nxloop
    CLOSE(10)

    ! Now we know how much memory we need for the temporary array
    ALLOCATE(exps(MAXL,MAXVAL(nxmax)))

    ! Open input file
    OPEN (10, FILE=filename, STATUS='old', ACTION='READ')    

    ! Skip first line
    READ (10,*) line

    nbasloop: DO 
       ! Read in line
       READ(10,'(A)', IOSTAT=sta) line
       IF (sta /= 0) THEN
          EXIT nbasloop 
       END IF

       ! Look for angular momentum shells
       DO l = 1, MAXL
          IF (str_lo_to_up(line(1:1)) == l_name(l)) THEN
             ! A contraction found.
             nbas(l,2) = nbas(l,2)+1
             ! Read in the number of exponents
             READ(line,*) dummystr, nxtmp, dummyreal

             ! Loop over the exponents
             DO i=1,nxtmp
                ! Read the exponent and the contraction coefficient
                READ(10,*) expn,dummyreal

                ! Check whether the exponent is already tabled
                found = .FALSE.
                DO j=1,nbas(l,1)
                   IF(expn==exps(l,j)) THEN
                      found = .TRUE.
                      EXIT
                   END IF
                END DO

                ! If it was not found, add it to the table
                IF(.NOT. found) THEN
                   nbas(l,1) = nbas(l,1) + 1
                   exps(l,nbas(l,1)) = expn
                END IF
             END DO
          END IF
       END DO
    END DO nbasloop

    CLOSE(10)

    DEALLOCATE(exps)

    IF(MAXVAL(nbas)==0) THEN
       CALL print_vertbar('!')
       CALL print_center('Error! Unable to find any basis functions in Gaussian''94 format!','!')
       CALL print_vertbar('!')
       CALL abort()
    END IF

  END SUBROUTINE read_gaussian_nbas
 !-----------------------------------------------------------

  ! Read in basis in Gaussian'94 format
  SUBROUTINE read_gaussian_coeffs(filename,X_l,C_l)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    ! Exponents and contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(OUT) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(OUT) :: C_l
    ! Helper arrays
    CHARACTER(LEN=MAXLEN) :: line, dummystr
    REAL(KIND=prec) :: dummyreal
    ! Read-in helpers for exponent and contraction coefficient
    REAL(KIND=prec) :: expn, contr
    ! Number of exponents and contractions
    INTEGER :: nx
    INTEGER, DIMENSION(MAXL) :: nxtabled
    INTEGER, DIMENSION(MAXL) :: nctabled
    INTEGER :: i, j, l
    ! Read-in status
    INTEGER :: sta
    ! Was exponent tabled already?
    LOGICAL :: found

    X_l=0.0_prec
    C_l=0.0_prec

    nxtabled=0
    nctabled=0

    ! Open input file
    OPEN (10, FILE=filename, STATUS='old', ACTION='READ')

    ! Skip first line
    READ (10,*) line

    nbasloop: DO 
       ! Read in line
       READ(10,'(A)', IOSTAT=sta) line
       IF (sta /= 0) THEN
          EXIT nbasloop 
       END IF

       ! Look for angular momentum shells
       DO l = 1, MAXL
          IF (str_lo_to_up(line(1:1)) == l_name(l)) THEN
             ! New contraction found.
             nctabled(l) = nctabled(l) + 1

             ! Read in the number of exponents
             READ(line,*) dummystr, nx, dummyreal

             ! Loop over the exponents
             DO i=1,nx
                ! Read the exponent and the contraction coefficient
                READ(10,*) expn,contr

                ! Check whether the exponent is already tabled
                found = .FALSE.
                DO j=1,nxtabled(l)
                   IF(expn==X_l(l,j)) THEN
                      ! Set contraction coefficient
                      C_l(l,j,nctabled(l)) = contr
                      found = .TRUE.
                      EXIT
                   END IF
                END DO

                ! If it was not found, add it to the table
                IF(.NOT. found) THEN
                   nxtabled(l) = nxtabled(l) + 1
                   X_l(l,nxtabled(l)) = expn
                   C_l(l,nxtabled(l),nctabled(l)) = contr
                END IF
             END DO
          END IF
       END DO
    END DO nbasloop

    CLOSE(10)

  END SUBROUTINE read_gaussian_coeffs

!-----------------------------------------------------------
  ! Write basis in Gaussian'94 format
  SUBROUTINE write_gaussian_basis(filename,X_l,C_l)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(IN) :: C_l

    INTEGER :: i, j
    INTEGER :: nx, nc
    INTEGER :: nprim

    INTEGER :: l

    nx=SIZE(C_l,2)
    nc=SIZE(C_l,3)

    IF(SIZE(C_l,2) /= SIZE(X_l,2)) THEN
       WRITE (*,*) 'Error in write_gaussian_basis: contraction coefficients do not correspond to exponents!'
       CALL abort()
    END IF

    ! Replace old file
    OPEN (20, FILE=filename, STATUS='replace')

    ! Write dummy element
    WRITE (20, '(A,I6)') 'E',0

    ! Loop over angular momentum values
    DO l = 1, MAXL
       ! Loop over contractions
       DO i = 1, nc
          ! Count number of primitives in this contraction
          nprim = 0
          DO j = 1, nx
             IF(C_l(l,j,i) /= 0) THEN
                nprim = nprim + 1
             END IF
          END DO

          ! Write out the shell
          IF(nprim>0) THEN
             WRITE(20,'(A,I5,F7.2)') l_name(l),nprim,1.0
             DO j = 1, nx
                IF(C_l(l,j,i) /= 0) THEN
                   WRITE(20,'(E20.12,E20.12)') X_l(l,j),C_l(l,j,i)
                END IF
             END DO
          END IF
       END DO
    END DO

    ! Write closing
    WRITE (20,'(A)') '****'

    CLOSE(20)

  END SUBROUTINE write_gaussian_basis

!-----------------------------------------------------------
  ! Write basis in Psi3 format
  SUBROUTINE write_psi_basis(filename,X_l,C_l)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: filename
    REAL(KIND=prec), DIMENSION(:,:), INTENT(IN) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(IN) :: C_l

    INTEGER :: i, j
    INTEGER :: nx, nc
    INTEGER :: nprim, np

    INTEGER :: l

    nx=SIZE(C_l,2)
    nc=SIZE(C_l,3)

    IF(SIZE(C_l,2) /= SIZE(X_l,2)) THEN
       WRITE(*,*) 'Error in write_psi_basis: contraction coefficients do not correspond to exponents!'
       CALL abort()
    END IF

    ! Replace old file
    OPEN (20, FILE=filename, STATUS='replace')

    ! Write dummy opening
    WRITE (20, '(A)') 'atomtype: "basisname" = ('

    ! Loop over angular momentum values
    DO l = 1, MAXL
       ! Loop over contractions
       DO i = 1, nc
          ! Count number of primitives in this contraction
          nprim = 0
          DO j = 1, nx
             IF(C_l(l,j,i) /= 0) THEN
                nprim = nprim + 1
             END IF
          END DO

          ! Number of primitives printed
          np=0
          ! Write out the shell
          IF(nprim>0) THEN
             ! Write out the type of the shell
             WRITE(20,'(A4,A1,A)', advance='no') " (",l_name(l)," "
             DO j = 1, nx
                IF(C_l(l,j,i) /= 0) THEN
                   ! Write out the exponent and the contraction coefficient
                   WRITE(20,'(A,E20.12,E20.12,A)',advance='no') "( ",X_l(l,j),C_l(l,j,i)," )"
                   np=np+1
                   ! If this was the last primitive, then close both brackets
                   IF(np==nprim) THEN
                      WRITE(20,'(A)') ")"
                   ELSE
                      WRITE(20,'(A)') ''
                   END IF
                END IF
             END DO
          END IF
       END DO
    END DO

    WRITE (20, '(A)') ')'
    

    CLOSE(20)

  END SUBROUTINE write_psi_basis

!-----------------------------------------------------------
END MODULE inout

